Excitable systems with noise and delay with applications to control: renewal theory 

approach 
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We present an approach for the analytical treatment of excitable systems with noise-induced 
dynamics in the presence of time delay. An excitable system is modeled as a bistable system with 
a time delay, while another delay enters as a control term taken after [Pyragas 1992] as a difference 
between the current system state and its state r time units before. This approach combines the 
elements of renewal theory to estimate the essential features of the resulting stochastic process as 
functions of the parameters of the controlling term. 
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I. INTRODUCTION 

Oscillations induced by noise are observed in many 
nonlinear systems with dissipation including neural net- 
works, lasers, chemical reactions, biomembranes, etc. 
(see [l| and references therein). At present, three differ- 
ent types of systems are being recognized that demon- 
strate noise- induced oscillations: bistable systems [H, 
systems close to Andronov-Hopf bifurcation [3[ , and ex- 
citable systems 

Much effort has been put into the development of the 
analytic description of excitable systems. The approach 
of the renewal theory [1] was successfully adopted to de- 
scribe the real excitable systems with continuous dynam- 
ics by a two-state @ , or by a three-state pi @ discon- 
tinuous stochastic process. In the frames of the renewal 
theory, it is assumed that the system is allowed to be only 
in a finite number of states, and the noise makes the sys- 
tem switch between the states in a random manner. The 
time the system spends in a given state before undergoing 
the next transition is called residence time in this state. 
For the validity of the renewal theory it is crucial that 
the distribution function of the residence time in a given 
state does not change in time. Therefore, in the case of 
an n-state system, its entire dynamics is determined by n 
distribution densities of the residence times, also known 
as residence time densities (RTDs). For the FitzHugh- 
Nagumo system the RTDs were explicitly calculated in 
using the Fokker-Planck equation approach suggested 
in 0]. Phenomenological three-state model of the ex- 
citable system with the arbitrary RTDs was discussed in 

Si- 

The situation becomes more complicated when a delay 
term is introduced into an excitable system, and the ran- 
dom process in it becomes non-Markovian. This problem 
arises for example in relation to the problem of control- 
ling noise-induced motion. 

Usually noise-induced oscillations possess a certain 
timescale that is dependent on the parameters of the ap- 
plied noise, e.g. its intensity. One possible way to de- 
fine this timescale for an excitable system is to estimate 



the power spectrum that would normally contain one or 
more peaks, and to take the inverse of the frequency of 
the highest peak to be the main period of oscillations. In 
[8| the idea was introduced to control the properties of 
oscillations induced merely by external noise by apply- 
ing time-delayed feedback force F(t) in the Pyragas form 
@: F(t) — k\x(t — t) — x(t)), where r is time delay and 
k is the feedback strength. It has been shown that the 
timescale and the coherence of noise-induced oscillations 
can be changed by adjusting solely the delay time in the 
controlling force. Moreover, an almost piecewise-linear 
dependence of the main period on the delay time was 
revealed. Remarkably, a similar behaviour of the main 
period was found in systems with timc-dclaycd feedback 
that were either excitable, or close to Andronov-Hopf bi- 
furcation [1, [l(| • For a van der Pol oscillator this phe- 
nomenon was explained by means of computation of the 
power spectrum of oscillations analytically: by using ei- 
ther linearization [Til . Il2l . [l3| , or mean- filed approxima- 
tion [H, [3. 

The theory of the excitable systems with time delay 
is still missing. Some progress was made in (l5j . where 
the impact of the delayed feedback on the current dy- 
namics of the system was considered in the mean field 
approximation. Here we propose an alternative approach 
based on the analysis of renewal processes with history- 
dependent RTDs. As mentioned above, the renewal the- 
ory can only be applied in the case when the RTDs do 
not change in time. This condition is obviously violated 
for the processes with time delay, in which RTDs are de- 
pendent on history. 

In order to overcome this problem, we introduce equi- 
librium RTDs by averaging over all possible histories. 
This approach is somewhat similar to the one suggested 
in [16[ for calculating the RTD in the case of stochastic 
resonance, where the noise-induced switchings between 
the two potential wells are considered in the presence of 
a weak periodic perturbation. Namely, in [161 ] the un- 
known RTD is obtained by averaging the known distri- 
bution of escape times (time it takes to escape from a 
certain well) over the known probability that the system 
spends a certain amount of time in the left well before en- 
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tering the right well. In contrast to this, we show that in 
the case of history-dependent renewal process, the equi- 
librium RTDs are given by the solution of an integral 
equation. This equation is derived for an arbitrary, and 
solved analytically for a moderate, delay time. 

In order to compare the analytic results of the mod- 
ified renewal theory with the numerical results for an 
excitable system, we design an excitable system from a 
bistable system with a non-symmetric potential and with 
a time delay, following the idea of the delay-induced ex- 
citability proposed in [l7|. To this system we add the 
controlling force F(t) as above. The parameters of the 
two-state model used for analytics are matched to the pa- 
rameters of the bistable system via the Kramers formula 
for the transition rates [18, 19, 20]. We employ the hybrid 
approach by combining the renewal theory results with 
the equilibrium RTDs of the history-dependent process in 
order to approximate the power spectrum of the noise- 
induced oscillations. Incoherence maximization due to 
delay is demonstrated for positive feedback strength k 
only. 



Switching 




X 



FIG. 1: Phase plane of FitzHugh-Nagumo system Eqs. ([T]) in 
an excitable regime. Dashed grey lines - null-clines, empty 
circle - fixed point, shaded area - the area from which an 
excursion can start. Stages of one full oscillation are shown: 
black dashed line - excursion, black solid line - refractory, 
dotted lines - switchings between the two branches (very fast). 



II. THE MODEL 

A. Excitability 

Before describing the way to construct an excitable 
system from a bistable system with delay following [TtI ] . 
we need to explain the concept of excitability. A popular 
toy model of an excitable system is a FitzHugh-Nagumo 
system 

/ x 3 \ 1 
X = [ X -T- y )-e> 
y = x + a + D((t), (1) 

where £(t) describes random fluctuations with Gaussian 
distribution, zero mean and unity variance. The null- 
clines of this system, which are the curves defined by 
x—0 and y—0 (assuming that D=0) are shown in Fig. [T] 
by grey dashed lines. Usually this system is considered at 
e<Cl in order to provide timescales separation described 
below. At a>l the system has a single stable fixed point 
(empty circle in Fig. [T]), and no oscillatory dynamics 
without noise (D=0). When Gaussian noise is applied 
(D>0), the behavior of the system changes drastically. 
Namely, while the value of D£(t) is small, the system os- 
cillates randomly around the fixed point during the wait- 
ing phase of duration T\y (note the horizontal plateaus 
of the realization of x{t) in Fig.[2ja)). But when the val- 
ues of D£(t) are large enough to throw the phase point 
into grey area in Fig. [TJ the system quickly tends to the 
right-hand branch of the cubic parabola along an almost 
horizontal path (lower dotted line in Fig. []}, since due 
to the smallness of e the horizontal component x of the 
phase velocity is much larger than its vertical component 
y. Then the system enters its excursion stage in which 



the phase point slowly crawls upwards along the right- 
hand branch of the parabola during Te time units (black 
dashed line in Fig. [T|): this motion is smeared by noise, 
but its velocity is almost unaffected by the latter, at least 
in average. As soon as the phase point reaches the top of 
the right-hand parabola branch, the vector flow swiftly 
carries it towards the left-hand branch (upper dotted line 
in Fig. [1]). Finally, the system enters its refractory stage 
as the phase point slowly crawls downwards towards the 
fixed point during time Tr (black solid line in Fig. [1]), 
again almost unaffected by noise in average. Then the 
process is repeated. 

A typical realization of a stochastic process x(t) occur- 
ring in system Eqs. (fTJ) looks similar to the profile shown 
in Fig. [DJa), where an x- variable from the FitzHugh- 
Nagumo system in the excitable regime is plotted. The 
cells in the lower panel show different stages of the pro- 
cess: white - waiting, shaded - excursion, patterned - 
refractory. One can single out two essential features of 
this motion induced merely by external noise. First, one 
can distinguish between very fast motion between the two 
branches of the cubic parabola, and slow motion along 
the branches. If e is very small, one can assume that the 
switching between the branches occurs instantly, i.e. the 
system can be only in one of the two states corresponding 
to the two branches of the parabola. Second, the dura- 
tions Te and Tjj of the excursion and of the refractory 
stages can be approximately regarded as independent of 
noise and constant in any event of a large excursion in the 
phase space, while the duration T\y of the waiting stage 
is completely determined by noise. The distribution den- 
sity of T\y depends both on noise and on the design of 
the system. 

The dynamics of the FitzHugh-Nagumo system in the 
excitable regime can be approximated by a two-state 
stochastic process s(t) which can take only two values, 
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Noisy bistable system with a single delay 
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FIG. 2: (a) Realization x(t) of an excitable system Eq. ([T]); (b) 
Stages of oscillations in (a) modeled by a two-state process; 
(c) Realization of a noisy bistable system with delay Eq. 
Lower panels of (a)-(c) show stages of the process: white - 
waiting with duration TV, shaded - excursion with duration 
Te, patterned - refractory with duration Tr. 



say, s = ±1 (Fig.[2Jb)). We assume that in the excursion 
stage the system is in state s = +1, and in the waiting 
stage in s = — 1. In addition, we assume that in the re- 
fractory stage s(t) does not change and assign s = — 1 
like in the waiting stage. Therefore, residence time in 
the state (—1) is equal to the sum of the refractory and 
waiting times, (Tr + TV). 

The dependence of the transition rate A on the param- 
eters of the FitzHugh-Nagumo system is rather complex 
and can be obtained only numerically. However, for a 
bistable system this information is readily available via 
the Kramers formula [H, [H, H(| • One can also neglect 
the intra-well dynamics of a bistable system and match 
the parameters of the latter with the parameters of a 
two-state system. 



B. Bistable excitable system with delayed feedback 

The idea of the present approach is to construct an 
excitable system from a bistable system by including a 
delay term into the latter, following the concept of delay- 
induced excitability proposed in pjj. Namely, in [l7| it 
was shown that a two-state system with a single time 
delayed feedback T can behave as an excitable one if the 
double-well potential describing the system is asymmet- 
ric. A characteristic feature of this model is that there is 
a locking of state, i.e. the transition from one state into 
another is allowed to occur not earlier than T seconds af- 
ter the moment of the previous transition. This feature 
is illustrated in Fig. [^c) where a realization of such a 
bistable system is shown, with the deeper potential well 
being located around the value of x = — 1. Here, one 



can recognize the same stages that occur in an excitable 
system (compare with (a)): waiting marked as a white 
cell in lower panel, excursion marked by shaded, and re- 
fractory marked by a patterned cell. Note, that in the 
bistable system with a single delay T the following condi- 
tion is automatically satisfied by the way of construction: 
Te = Tr = T . As a result, the distribution densities of 
the waiting times in both states are shifted to the right 
by T as compared to the case without delay. 

We would like to assess the effect of the controlling 
term F(t) on this system, in which the second delay ap- 
pears. The second delay term is due to the delayed feed- 
back force. 

The dynamics of the bistable system with two time 
delays is described by the equation 



x(t) 



dU[x, xt] 
dx 



(2) 



where xt and x T denote the retarded variables x(t — T) 
and x(t — t), respectively; T is the fixed excursion and 
refractory time, r is the delay time of the controlling 
force, k is the strength of the controlling force, is a 
Gaussian noise, and D is the noise strength. 

The potential U[x, xt] in Eq. ([2]) is chosen to have two 
wells with minima located at x = ±1 that are separated 
by a barrier with the maximum at x = —0.3, so that 

4 2 

U[x,x T ] = — - y-(0.1 + 0.2a;T)(a; 3 -3ar), (3) 



dU[x, xt] 
dx 



(x 2 -1)(x-0.3-0.6xt). 



Following [l8( , we assume the overdamped case when we 
can neglect the intra-well dynamics. It is well known 
that this approximation is valid for noise intensities much 
smaller than the height of the potential barrier, i.e. for 
D < Af/cff. 

The special choice of the position of two minima does 
not affect the results presented below. 

The evolution equation Eq. ^ can also be rewritten 
in terms of an effective potential U c s which includes the 
controlling force 

x 2 

U c s[x, xt, x t ] = U[x, xt] + k— — kx T x. (4) 

We will be considering small values of feedback strength 
k for which xt and x T can be either in the left or in the 
right well. With this, the values they are taking are close 
to —1 or to +1, respectively, therefore we single out two 
states of the system that we denote as (—1) and (+1). 
In the analytical calculations below we will substitute 
the exact instantaneous values of xt and x T by their 
approximate values ±1. 

In Fig. [3] the effective potential U e s is shown for fixed 
xt and x T that take values ±1. Consider all stages of 
one oscillatory cycle. Start from xt ~ +1- From Fig. [3] 
it is clear that at any x T there exists only one well (left) 
and the particle trapped in it will remain there until xt 
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FIG. 3: Effective potential U e g{x,XT, x T ] Q for fc = 0.05 and 
fixed xt and x T as in the legend. 



switches to approximately (—1). This is the refractory 
phase with duration Tr. At xt~— 1 the right well ap- 
pears and the waiting phase begins, during which it is 
possible to jump from the left to the right well. Transi- 
tion rate A from (—1) to (+1) becomes dependent on x T : 
for k > it is larger when x T ~ +1, i.e. when the left 
well is shallower. After the particle jumps into the right 
well, the system is in the excursion phase. Note that at 
xt ~ — 1 the right minimum is much deeper than the 
left one. Therefore, the particle trapped in it will remain 
there until xt changes from (—1) to (+1), i.e. until the 
right well vanishes. After that the particle jumps from 
state (+1) to the state (—1) and the cycle repeats again. 
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TABLE I: Transition rates A between the states of the 
bistable system with delay at different values of xt and x T - 



The information on the transition rates is summarized in 
Table [U Here we denote by (p + q), q> —p, the transition 
rate from the state (—1) to the state (+1) if x T w +1 and 
xt ~ — 1 and by p, p > 0, the the transition rate from the 
state (—1) to the state (+1) if x T w —1 and xt ~ — 1. The 
transition rates are calculated from U r according to the 
Kramers theory [l8[ as follows 



« \J d xx U c ff(x m , 1, l)d xx U c f[(xo, -1, -1) 

Z7T 



exp 



D 



p + q = —\/-d xx U e s(x' m ,-l,l)d xx U eS (x' ,-l,l) 

Z7T V 



exp 



AC4 
D 



(5) 



where xq and x m are the positions of the maxima and the 
left minima of the effective potential U c g at xt = x T = 
— 1 and x' and x' m are the same quantities at xt = — 1, 
x T = + 1. At/off and AU^ S are the potential differences 
between the maxima and the minima. Note, that due to 
presence of the feedback term k(x T — x), the values of 
x m and x' m are not exactly equal to —1, and xo is not 
exactly —0.3. However, in analytic calculations we sub- 
stitute x — x'q, x m and x' m by their approximate values 
—0.3, —1 and —1, respectively. Throughout the paper we 
compare the analytic results derived from the two-state 
model with the results of simulation of a bistable system 
with two time delays. 



III. TWO-STATE MODEL: 
RESULTS 



ANALYTIC 



One should notice that the two-state model is not de- 
rived from the continuous bistable system Eq. ([2]) with 
linear feedback term. It can describe a larger class of 
systems with similar properties, e.g. systems with non- 
linear feedback term, provided that the transition rates 
depend on the delay time according to the phcnomeno- 
logical rule provided by Table (jj) . 

The continuous random process x(t) is approximated 
by the discrete random process s(i) = ±1 with infinitely 
fast (discontinuous) transitions from one state to the 
other. In order to match the parameters of the two-state 
model with those of the bistable system Eq. @, the ex- 
cursion and the refractory times must be set equal, i.e. 
Te = Tr pjj • However for the sake of generality the an- 
alytic results presented below were obtained for the case 
of Te 7^ Tr. In what follows by £ we denote the time of 
residence in a certain state: from the context it will be 
clear what state or phase is referred to. 

Regarding the dependence of the transition rates on 
the history (Table [£]) we conclude that the RTD ip + (£) 
in the state s = +1 is given by the delta- function 
tp+{£.) = <5(£ — Te). This assumption is justified if there 
is a strong time scale separation meaning that the transi- 
tion from the excited to the non-excited state and back- 
wards occurs almost instantaneously. Moreover, the RTD 
in the state s — — 1 is zero for the first Tr seconds after 
the transition from the state s = +1 to the state s = — 1. 
Consequently, the model contains only one unknown ob- 
ject, namely the distribution density tp-(0 of the waiting 
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FIG. 4: (a) Sketch of the profile of the two-state process s(t) = 
±1 on the interval of time [to - T E - 2T R ;t ]. (b) Three 
possible forms of the survival probability P_ (£) in the waiting 
phase for r < Te + 2Tr. Here, £ is the duration of the waiting 
phase. P- (£) is shown for three different values n of the delay 
time t chosen arbitrarily from the intervals indicated in the 
legend. Filled circles on the £-axis indicate points where the 
change of the slope of the corresponding P-(f) occurs. 



time. It is related to the RTD ^_ (£) in the state s = — 1, 
via 



MO 



0, £e[0;T fl ] 
^-(e-Tfl), [T R ;oo) 



(6) 



A. Small delay times r: non-variable history 

We start with the case of small delay times from the 
interval r G [0;T E + 2T R ]. The probability P_(£) of 
survival in the state (—1) during the waiting phase is 
given by the solution of the equation 



M0 = 



gMO 



= -A_(0M0, 



(7) 



where A_(0 stands for the transition rate from s = — 1 
to s = +1. If A=const, which is valid for the bistable 
systems without delays, the statistics of the waiting times 
is exponential [3] with the distribution density ip given 
by 



^)=Aexp(-A£), 



(8) 



where the constant transition rate A depends on the sys- 
tem parameters. The mean waiting time (TV) is then 
given by (TV) — 1/A. However, due to delay there ap- 
pears a discontinuity in the transition rate A_(£) (see 
Table 1), therefore the solution of the Eq. ([7]) on the in- 
tervals where A_ (£) is constant must be normalized in 
such a way that the survival probability P_ (£) remains 
continuous. 

Suppose that at time t — t (Fig. UK a)) the refractory 
phase has just finished and the system is in the very 
beginning of the waiting phase, implying that the waiting 
duration is equal to zero, £ = 0. Fig. |D(a) shows the 
profile of the two-state process on the interval of time 
t £ [to — Te — 2T R ; to}. It is clear that no other profile 



s(t) is possible on this interval of time. This means that 
despite the fact that the solution of Eq. ([7]) depends on 
r, this dependence remains tha same during any waiting 
phase of any cycle, as long as r is less than Te + 2T R . 
Hence, the distribution density of the residence times in 
the state s = — 1 does not change in time and the renewal 
theory can be applied. Three different cases should 

be considered separately for r <Te + 2T R . 

Case 1: r G [0;Tr]. The solution of Eq. (7J is then 
given by 



P_(0 = e~»* £e[0;oo). 



(9) 



This solution is shown by the dotted line in Fig. Hfb). 

Case 2: r G \Tr\Te + T R ]. In this case the survival 
probability reads 



MO 



,-(p+q)t 



£e [0; 



e^e-iir-Tn)^ t e [ T -T R ; R oo). (10) 



This solution is shown by the solid line in Fig.[3Jb). 

Case 3: r £ [T E + T R ; T E + 2T R ] . 



MO 



(11) 

1, Ce[0;r-T B -T fl ] 
e -?e e ?(r-T B -T R ) ) £e[T-T E - T R ; 

r-T R ] 

e~ qTE , ee[r-T B ;oo). 

Dashed line in Fig.[4jb) corresponds to P_(0 in the Case 
3. 

We can define the main period T ma in as r ma in = 
27r/fi max , where the frequency f2 ma x corresponds to the 
absolute maximum of the power spectrum of oscillations. 

In [8[ it has been shown numerically for FitzHugh- 
Nagumo system that T ma i n depends almost piecewise- 
lincarly on r. It was assumed there that the piece- wise 
linearity in the main period vs delay time is a univer- 
sal phenomenon that should occur in any system driven 
by noise including non-excitable systems near Andronov- 
Hopf bifurcation. Here, a similar dependence of T ma j n on 
r will be demonstrated for system Eq. © analytically 
already for r G [0;T S + 2T R ] and fc<0. From Eqs. ©- 
(fTTj) one can calculate the power spectrum of oscillations 
using the result of the renewal theory 0, HH 



S(u) 



2(Ax) 2 1 Re I [l-tt__M][l-sM»q,)] ^ a2) 



(T+) + (T_) u>* 



1 - ^-(iu>)tl) + (iu)) 



where ip+(iu))=exp (— iujTe) is the Laplace transform of 
ip+, *f>-(iuj) is the Laplace transforms of , F_ from Eq. ((6]), 
Ax=2 and (T±) are the mean residence times in the 
states s—+l and s=— 1 of the two-state process, re- 
spectively. It is clear that (T + )=Te and (T-)=Tr + 

Jo°°^-(OC 

We compare the analytic power spectrum (|12p with the 

numerical power spectrum of the original bistable sys- 
tem @ calculated with the parameters k=— 0.05, T=50, 
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1 L5 ~~ 2 2.5 3 

x/r E 

FIG. 5: (a) Numerical power spectrum of the bistable system 
for T — 50 and three different r close to r c . (b) Analytic 
power spectrum (112[) for the same parameters as in (a) and 
Te=Tr— 55. (c) Numerical (circles) and analytic (line) main 
period of the oscillations. 



D=0.053. Note that any switching event in the bistable 
system has a finite duration, i.e. the change of the co- 
ordinate x from (+1) to (—1) takes certain time. This 
leads to the increase of the duration of locking. There- 
fore to successfully match the parameters of the bistable 
system @ to the parameters of the two-state model we 
set Te=Tji=T+A, where A is the time needed for the 
particle to descend from the maximum of the effective 
potential U c s into its minimum. A is approximately set 
to 5 for |fc|<0.05. 

For the delay time r between 2Te and 3Te the power 
spectrum has two maxima as shown in Fig. EJa). For r 
close to 2Te the left maximum is higher than the right 
one (dashed line in Fig.[5fa)). As r increases from 2Te to 
3Te, the left maximum of S(oj) decreases and the right 
maximum increases. At critical t c the two maxima have 
equal heights as shown by the solid line in Fig. 03a) . For 
t>t c the right maximum is higher than the left one as 
shown by the dotted line in Fig. Ufa). Fig. H^b) shows 
the analytic power spectrum for the same parameters as 
in Fig. GJa). 

The main period increases almost linearly with r for 
r<r c , at t—t c the main period drops discontinuously and 
for t>t c it increases almost linearly again. This is shown 
in Fig.[5fc), where numerical results (circles) are com- 
pared with main period computed from the analytical 
power spectrum ([12"]) (solid line). 



B. Large delay times r: Equilibrium RTDs. 

Now consider the case of the delay times larger than 
(Tb+2Tr). This situation is qualitatively different from 



the one discussed in Section [ill Al because for large delay 
times the history is variable. In other words, there are 
many possible profiles of the two-state system on the in- 
terval of time [to — t;4q], where to is the time moment 
when the current refractory phase has just finished. This 
situation is sketched in Fig. [SI where the time moment to 
is indicated by the filled circle. The distribution of the 
residence times in the state s=— 1 changes in time and 
now the results of the renewal theory cannot be applied 
directly. To overcome this problem we introduce the con- 
cept of equilibrium RTDs in the sense of the averaging 
over all possible histories. The equilibrium RTDs can be 
computed from the known instantaneous RTDs as shown 
below. 

Consider a certain history s(t) of the two-state stochas- 
tic process s(t)=±l on the interval of time t £ [to — t; to] 
inside which there are N T pulses. Note, that inside an 
interval of length r there can be not more than N T 
pulses, where N t =[t /{Te + Tr)] and [•■•] denotes an 
integer part of the number. Number these pulses with 
index j changing from 1 to N T < N T , with j=l cor- 
responding to the pulse closest to to- Denote by Uj 
the duration of the waiting phase that precedes the j- 
th pulse. If we know the number N T of pulses inside 
t € [to ~ T ! to] together with the ordered durations of 
the waiting phases Uj, j = 1,2, ...,N T , we can un- 
ambiguously reconstruct s(t). Let us consider all pos- 
sible histories §k(t) and form an iV r -dimensional "his- 
tory" space with vectors Ufc=(ttfe i x, Wfc,2j ■ ■ • , Uk,N T )- 
If the number of pulses N Ti k for the given history §k(t) 
is less than the largest possible number N T , the redun- 
dant coordinates Ufcj, j—N T ^+l, ■ • • , N T are set to zero. 
Namely, if there is only one pulse on the given interval, 
we set U).,i=(t—Te-T r ) and Uk,j= 0, j=2,...,N T . If 
there are two pulses, we set Uk,2=(T— 2Te~ 2Tr— Uk,i) 
and Ukj—0, . . . , N T , etc. Therefore, the durations 

Uk,j cannot be larger than (r — Te — Tr) by the way of 
construction. Obviously, is in one-to-one correspon- 
dence with the set of all possible histories Sfc(t). 

Next, split the space into M equal A^-dimensional 
cells Ii, I — 1, 2, M of the form Ij=[ttj ) i ± rj/2] x [tij )2 ± 
T}/2] x ... x [tt|r j jv. r ±?7/2], where t] is the side length of each 
cell. can be represented as a unity of A^-dimcnsional 
cells as follows: 

HI = {hUl 2 U...Ul M }. 

Consider one cell Ik- One can assume that all histo- 
ries corresponding to the points inside this cell are ap- 
proximately the same. Therefore, the waiting times that 
start immediately after to and whose histories are ap- 
proximately Ufc , have approximately the same RTDs ip^ k 
which are solutions of Eq. ([7]) . 

Now consider some sufficiently long realization of the 
two-state stochastic process s(t). Suppose we are inter- 
ested in some function / of the randomly changing wait- 
ing time £, and we wish to compute an average of /. 
There are two kinds of averages: over time and over the 
ensemble of realizations, which coincide if the underlying 
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FIG. 6: Schematic representation of the configuration of 
pulses in the case of r G [T E + 2Tr; 2T e + 3Tr]. (a) The 
durations of the two latest waiting times are Tw = 0. (b) 
Durations of the latest waiting time is Tw = u > 0. (c) 
Durations of the latest waiting time is Tw > u. 



process is ergodic. The time average (/) can be calcu- 
lated from a single realization as follows 



1 



JV 



N^ca N L — ' 
t=l 



(13) 



where N stands for the number of pulses in the realiza- 
tion. In a sufficiently long realization one can find a large 
number rik of histories from the cell Ik ■ The waiting times 
occurring after these histories have (approximately) the 
same distribution ?/;" fe . In an infinitely long realization 
when N— >oo, one can find all histories from the same 
cell, i.e. rife — » oo, with every history being found an in- 
finite number of times. We now regroup the summands 
in Eq. (TIB")) by collecting inside each bracket number k 
the values of / of the waiting times £ occurring after the 
history from the same cell centered at u^.. We denote the 
respective values of / as /"* , I = 1, . . . , n^. Also, each 
bracket number k is divided and multiplied by n^. 

</> = ^fe^ J/" 1 +/" 1 + --- + g] 



n 2 1 



N n 2 <. 



[/r+/ ; 



u 2 

2 



-CI 



divided by the total number N of pulses in any given re- 
alization of the stochastic process, in the limit as N — > oo 



C q (0^= lim §. 



(15) 



Clearly, the factors n^/N in (fT4"]) are the probabilities of 
the history to be in the cell Ik- Denote these probabil- 
ities by V(uk) (drj) NT , where V(u) is the corresponding 
distribution density. From the general considerations it 
is clear that V(vl) depends solely on ?/>° q . The particular 
dependence of V(uk) on ■ip cq is to be determined sepa- 
rately for any given history 

The terms ^ [/"" +f% k +... + f^] in Q3J) are the 



time averages of /(C), for £ with the history determined 
by the cell centered around Ufc and with distribution den- 
sity ip- k ■ Now, since the process is assumed to be ergodic, 
a time average is equal to the ensemble average. The lat- 
ter can be calculated with the knowledge of the RTD ■0" fc 

for the given history as J °° V'- fc (£)/(£) d£- 

Finally, in the limit M — > oo and t] — ► 0, the equation 
(11411 can be rewritten in the form 



(/> = / iT(0/(£K 

Jo 

= f V{n)du 1 du 2 ...du NT n^{§f{0de. 
JueHi Jo 



(16) 



Regroup terms in Eq. ([16]) as follows 



o= / m<% -«)+/ v(u)i>»(0du 1 du 2 ...du^ 

10 I JuEHT, 



Note that Eq. (flu]) holds for an arbitrary integrable func- 
tion /, which means that the expression in brackets is 
equal to zero. We therefore arrive at the following inte- 
gral equation for the unknown equilibrium RTD 



C q (0= / V(u)^)d Ul du 2 ...du NT . 



(17) 



Equation (fT?)) allows us to treat the renewal processes 
with non-identically distributed waiting times. It gener- 
alizes the results of the renewal theory [4J . 



riM 1 
N n M - 



C. Equilibrium RTD for r G [T E + 2T R ; 2T E + 2T R ] 



Introduce the equilibrium RTD V , - <1 (£) in the sense that 
"0- q (O^ gives the probability for the waiting time to 
have the duration in the interval [£; £ + d£\. The equilib- 
rium RTD does not depend on history, i.e. ^l q (£) dt; is 
understood as the number N$ of the waiting times with 
the duration Tw that fall within the interval [£; £ + d£] 



We now use the derived equation (fT?)) to calculate the 
equilibrium RTD ip cq in the case of r G [T E + 2T R ; 2T E + 
2T R ]. The reason for choosing this interval for r in- 
stead of the whole interval [T E + 2T R ; 2T E + 3Tr], where 
the history space is one-dimensional is the following. 
If r G [T E + 2T R ; 2T E + 2T R ] the history can contain 
only one complete pulse and one half-complete pulse (see 
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Fig. [6J, depending on the duration of the previous wait- 
ing phase Tyf . Therefore, this situation is simpler than 
the case of r G [T E + 2T R ; 2T E + 3T R ], where the history 
can contain up to two complete pulses. Denote the dura- 
tion of the previous waiting phase by u and consider two 
possible cases. 

Case 1: For u G [t— Te — 2T R ; oo) the history contains 
one complete pulse. In this case the solution of Eq. (JT]) 
is given by 

P°(0 







(18) 


f 1 ' 




- T B - T fl ] 


1 e -qi+q(r~T E -T R ) ^ 


Ce [r- 






T — 


Tn] 




[r- 


T R ;oo). 



Case 2: For u G [0; t — Tg — 2Tr] the history contains 
one complete pulse and one half-complete pulse. In this 
case the solution of Eq. ([7]) is 



(19) 



x < 



e" 


5 


[0; 










T — 


T B 


- 2T fl - u 


e" 


-g(r-T fi -2TH-^) 


£e [t- 


T B 


-2T R -u 






T — 


T B 


-T R ] 


e" 




fe [r- 


T B 


-T R ; 






T — 


3«] 




e" 


-g(T-2T R -u) 






oo). 



The connection between the probability density V (it) and 
■0l q is straightforward 



V{u) = O). 



(20) 



Because u is constrained to the interval [0; r — Te — 2T R ] , 
the relation Eq. (|17p is an equation for the unknown 
"0- q (f ) on ly on the same interval of £ G [0; r — Te — 2T R ]. 
For £ > (t — Te — 2T r ), Eq. (TIT)) is no longer an equation, 
as will be shown below. 

For £ G [0; r — Te — 2Tr] we divide the integration 
interval {/ according to the restrictions on £ and it as in 
and rewrite Eq. (JT7J) as follows 



1 - 



— T B -2T n 



-(p + q)e- ip+,l) t 



ip e _?(u)duj (21) 

t-T e -2Tr-£ 





-T E 



-2T r 



%l)°3(u) du 
e qu ^ c 3{u) du. 



It-T e -2T r -^ 

We look for the solution of Eq. (j2"Tj) in the form 

•0l q (t) = Ae-^ +q ^, (22) 



where A is some unknown constant. Plugging (122[) into 
(|2"TT) and comparing the coefficients of the exponents e~ p ^ 
and e~( p+9 ^ yields the constant A 



A = 



pip + q) 



p 



+ qe -(p+q)(T 



-2T R ) ' 



(23) 



On the interval £ G [r - T E - 2T R ; r - T E - T R ] Eq. |(T7|) 
becomes 



1 - 



t-T e -2Tjj 



-pe" p5 e _<?(T_Ti5_2TR) 



T-T E -2T n 



e"Vl q (24) 



which is no longer an equation, because the r.h.s. of 
Eq. depends on the known i/>l q on the interval £ G 
[0; t-T e - 2T R ] . Similarly, the equilibrium RTD V_ on 
the interval £ G [t — Tg — Tr; r — T R ] is determined as 



pt—Te — 2Tr N 

i/^ q (u) du 



(25) 



1 - 



+ (p + g)e-( p+IJ )«e 9T « 



r-T E -2T R 



e 9 "V! q (it) du. 



Finally, the equilibrium RTD i/?l q on the interval £ G 
[r — Tr; 00) reads 



pe -p? e -?(r-2T n ) 



t-T e -2T r 



ip_(u) du 





t-T e -2Tr 



e 9U Vl q (u) du. 

(26) 



Taking the integrals in Eqs. (|23 )) -(|2l)] ) with Vl q on the 
interval £ G [0;t-T e - 2T r ] from Eq. (J22J), one obtains 
the equilibrium RTD on the whole interval 



«) = A 
( e -(p+q)t 

q(T-T E -2T R )-pt 



X < 



(P+q) c -(p+q)£ c qT„ 
e _g( T _ 2T R )-p£ 



(27) 

£G [0;t-T e -2T r ] 
£ G [r - T E - 2T R ; 

t -T E - T R \ 
£G [t-Te-T r -t-T r ] 
£ G [t - T R ;oo). 



The knowledge of the RTD allows us to compute the av- 
erage of any given function g(s) of the stochastic variable 
s(t) 



(<?> = 



g (-l)(T_) + g(+l)T E 



(T- 



(28) 



where (T_) = T R + / °° uiM q (ii) du for r G [Te + 
2T R ;2T E + 2T H ]. For the values of the delay time 
t < Tg + 2Tr the equilibrium RTD must be replaced by 
the corresponding i/>_ (u) from Eq. {2} with P_ (it) from 
Eqs. CD) or (nj. 

We present here the results of the comparison of the an- 
alytic formula Eq. (j2"5)) with the simulation of the original 
bistable system @. Figs. Efa)-(b) show the average (x) 
vs delay time r for D = 0.043, jfe = 0.05 and T> = 0.053, 
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FIG. 7: Comparison of the average x vs delay time calcu- 
lated numerically (circles) with the average (s) of the two- 
state process (|28|) (lines). Parameters are (a): D = 0.043, 
T = 50, T E = T R = 55, k = 0.05 and (b): D = 0.053, T = 50, 
T e =Tr = 55, fc = -0.01. 



k = —0.01, respectively. Analytic result (j2"5j) (lines) is 
compared with (x) calculated numerically (circles). As 
in Fig. [5] the durations of the excited and the refractory 
phases are Te = Tr = 55. 

For positive feedback strength (Fig. [7(a)) the average 
(x) increases with the delay until r = 2Te then it de- 
creases when r is in the interval r 6 [2Te; 3Tb] and after 
that it increases again when r G [3Te;4Te]. The be- 
haviour of (x) on r is exactly the opposite for negative 
feedback strength k = —0.01 as it is shown in Fig. [7(b). 
We do not plot the variance (x 2 ) because it is indepen- 
dent of delay: according to Q28p. for the stochastic pro- 
cess s = ±1 the variance is constant ((s 2 ) = 1). 

Note that the non-monotonous dependence of the av- 
erage (x) on the delay time is qualitatively similar to 
the dependence of the variance (x 2 ) on delay time which 
was derived in [22j for the van der Pol oscillator near the 
Andronov-Hopf bifurcation. This similarity shows again 
that the properties of the noise- induced oscillations in the 
excitable systems are closely related to the properties of 
non-excitable noise-driven systems near bifurcations [|[ . 



IV. POWER SPECTRUM IN THE MEAN FIELD 
APPROXIMATION 

The knowledge of the equilibrium RTD does not allow 
us to compute the power spectrum of noise-induced oscil- 
lations. According to Q the power spectrum is given by 
the Fourier transform of the average product of proba- 
bility currents (j{t)j(t')}, where the current j(t) is a sum 
of 5-like pulses occurring at the moments t n of switch- 
ing: j — 2^2 n (— l) n 5 (t — t n ) . Any switching time t n 
is the time moment, when the transition from the state 
s = — 1 to the state s = +l or backwards occurs. The 
amplitude of the current j is given by +2 in case when 
s = — 1 changes to s = +1 and by —2 otherwise. Unfortu- 
nately, the average (j(t)j(t')) can not be represented in 
the form (fT4"]) and therefore the correlation function can 
not be written in terms of the equilibrium ip 1 ^ 1 . 

One possible approximation for the power spectrum 
is an analog of the mean field approximation, when the 
known result for the power spectrum from the renewal 




FIG. 8: (a) Circles represent the main period Tmain of the 
oscillations in units of Te = Tr = 55 calculated numerically, 
solid line shows the corresponding analytic counterpart com- 
puted from (|12|) and (|27l) . (b) Zoom of the region in (a) near 
r = 3Te- (c) Numerical power spectrum for r near the crit- 
ical t c . Dashed line shows power spectrum at r = r c . (d) 
Analytic power spectrum for Te = Tr — 49. The inset in (d) 
shows zoom of the area around the maximum of S(u). 



theory [1, [2l| is calculated with the equilibrium -0l q ob- 
tained above. 

We compare this version of the mean field approach 
with the simulation of the bistable system @ in Fig. [5] 
For D = 0.043, k = 0.05 and the rest of the parameters 
as in Fig. [JJ in Figs. H(a,b) the main period T ma i n of 
the noise-induced oscillations is given as a function of 
delay time r. Both T ma j n and r are shown in the units of 
Te = Tr. Circles in Figs.[8(a,b) correspond to the results 
of simulation, lines show the main period computed from 
the analytic expression for the power spectrum (|12p with 
the equilibrium RTD (|27p . For simplicity we show only 
the part of the analytical curve for r £ [3Tg; 4Tg]. As we 
see, the mean-field analytics predicts correctly the critical 
delay time where the first branch- switching occurs. 

In Fig. [Sic) the numerical power spectrum is shown 
for three different values of r close to its critical value 
T c- Fig. [8(d) shows the analytic power spectrum for the 
same parameters as in Fig. [5(c) and Te = Tr = 49. 



V. INCOHERENCE MAXIMIZATION DUE TO 
DELAY 

A. Bistable excitable system 

The RTDs ^_ derived in Section IIII1 along with the 
equilibrium RTD (|27p , allow us to compute the coefficient 
of variation R [l|] , which serves as a measure of coherence 
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FIG. 9: (a) Coefficient of variation R vs noise strength D for 
k = 0.05 and Te = Tr = 55. (b) The same as in (a) for the 
negative feedback strength k — —0.05. 



of spiking 



R 



y/(Tl) (T-)< 

T E + (T_) 



(29) 



The smaller is the coefficient of variation R the more 
regular is the spiking. The effect of coherence resonance 
[23j | manifests itself in the appearence of a minimum in 
the dependence of the coefficient of variation on the noise 
strength D. The opposite effect to coherence resonance 
is regarded as incoherence maximization [l[. It appears 
when there is a local maximum in the dependence of R 
on D. 

Here we show that in the bistable system presented 
in the SectionHTl delay induced incoherence maximiza- 
tion can be observed for positive feedback strength k. 
In terms of the two-state model with history dependent 
transition rates Tab. (p}, positive feedback strength k cor- 
responds to the enequality p + q > p, implying that 
the probability of transition from the state s = — 1 to 
the state s = +l is larger if s T = +l. Since the dura- 
tions of the excited Te and the refractory Tr states are 
fixed, the dependence on noise is realized through the 
duration of the waiting phase T\y only. Obviously, in 
the limit of large noise strength (D — »oo), the duration 
Tw becomes vanishingly small, leading to almost reg- 
ular spiking with the mean interspike interval given by 
Tb + Tr. For vanishing noise (D = 0) the transition prob- 
ability is neglidgibly small resulting again in the regular 
spike train. Therefore, for some finite noise strength we 
can expect to observe a local maximum of the spike in- 
coherence. 

Using the Kramers relation (JSJ [l^, H(| between the 
transition rates p and q and the noise strength D, we plot 
in Fig. © the coefficient of variation R as function of D. 



In Fig. [9ja) R is shown for positive feedback strength 
k = 0.05, Te = Tr = 55 and different delay times t as 
indicated in the legend. With no feedback (r < Te) the 
coefficient of variation decreases monotonically with D. 
However, at any t larger than, and close to, Te a maxi- 
mum appears at a certain noise strength D. This is an ev- 
idence of incoherence maximization induced by the delay. 
If the feedback strength is negative k = —0.05 (Fig.^b)) 
the incoherence maximization is absent at least for the 
delay times t < 4Tg. 



B. Comparison with the FitzHugh-Nagumo system 

To demonstrate incoherence maximization predicted in 
the Section (|V A[) , we use the FitzHugh-Nagumo system 
with non-linear delayed feedback introduced through the 
activator x into the equation for inhibitor y 



II 



K[x T 



oat), 



(30) 



where x T = x(t — r), e = 0.01 and a = 1.1. 

The form of the feedback term in Eqs. (f3"0"l) was chosen 
in such a way, that the noise-induced dynamics in the 
waiting phase is similar to that of the two-state model 
introduced in the Section (|III|1 . 

It should be emphasized that the FitzHugh-Nagumo 
system with linear feedback cannot be approximated by 
a two-state model with the transition rates as in Ta- 
ble (|T|. To see why this is so, consider a typical pulse 
train given by variable x as function of time. This is 
shown in Fig.fJJa). Assume that the system is currently 
in the waiting phase, x(t) £ waiting phase, and assume 
further that r seconds ago the system was in the re- 
fractory phase, x(t — r) e refractory phase. Since the 
value of the state variable in the waiting phase is differ- 
ent from that in the refractory phase, unlike in the two- 
state model, we conclude that the feedback term which is 
given by k[x(t) — x(t — r)] is not zero. However, this term 
becomes zero (up to fluctuations whose order is given by 
the noise strength D) if r seconds ago the system was in 
the waiting phase. 

This means that the current transition rate A from the 
waiting phase to the excited phase changes depending on 
whether r seconds ago the system was in the waiting or 
in the refractory phase. This contradicts the assump- 
tion that the transition rate depends on r according to 
Table. ©. 

The comparison between the FitzHugh-Nagumo model 
and the two-state model is possible only in the case of 
nonlinear feedback, e.g. like in Eqs. (|30p and strong time 
scale separation, i.e. e <C 1, when the concept of the tran- 
sition rates can be applied. The nonlinear feedback term 
ensures that during the waiting phase (x(t) ps —1.1) the 
value of [x r — x] 2 becomes significant only if r second ago 
the system was in the excited phase, i.e. if x T ~ 2. On 
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FIG. 10: (a) Coefficient of variation R vs noise strength D 
for the FitzHugh-Nagumo system Eqs. (|30[) with fixed delay 
time t — 4 and the feedback strength K as indicated in the 
legend, (b) Coefficient of variation R vs noise strength D for 
the two-state model with the parameters Tr = 2.5, r = 4, 
and different feedback strength k given in the legend. 



the other hand, if x T belongs to the refractory phase, the 
term [x T — x] 2 is negligibly small, so that the transition 
rate from the non-excited to the excited state is the same 
as in the original system without the feedback. Conse- 
quently, the transition rate is modified by the feedback 
only if x T belongs to the excited state. 

It is easy to see from Eqs. ([50]) that negative values of 
the feedback strength K effectively decrease the transi- 
tion rate, whereas the positive values of K increase the 
probability of transition. Therefore, based on the predic- 
tions made in the Section (jV A[) , we conclude that delay- 
induced incoherence maximization should be observed for 
positive K. 

This result is confirmed by numerically computing the 
coefficient of variation R vs noise strength D for the 
FitzHugh-Nagumo system Eqs. (|30|) at fixed delay time 
r = 4, as shown in Fig. (fT0")) fa). For negative K and 
K = 0, R decreases monotonically with D, however for 
positive K the coefficient of variation reaches a local 
maximum confirming incoherence maximization at cer- 
tain noise strength D. 

To match the parameters of the two-state model with 
those of the FitzHugh-Nagumo system Eqs. (f30|) . we set 
Tr = 2.5, Te = 0.5, t = 4 and plot R vs noise strength D 
in Fig. (fH))) (b) at different values of the feedback strength 
k (see legend). We see that the behaviour of the coeffi- 
cient of variation R on D in the two-state model quali- 
tatively coincides with the numerical results obtained for 
the FitzHugh-Nagumo system Fig. ([TO]) (a). 



VI. CONCLUSION 

To conclude, we presented a two-state model of an ex- 
citable system with time-delayed feedback. In this model 
the state variable s takes only two values s = ±1 and the 
transition probability from one state into the other de- 
pends on the history of the process in a given way. To 
compare the results derived for the two-state model with 



the properties of a real excitable system we consider a 
bistable system ([2]) with the effective potential (Q| which 
contains two delay times. One of them is fixed and is used 
to model the excitability and the other one is assigned to 
the delay time of the controlling feedback force. 

Assuming that the durations of the excited phase and 
the refractory phase are noise- independent, we conclude 
that the only unknown and in general history-dependent 
quantity is the residence time density (RTD) ip-(0 °f 
the waiting phase. We show that for the delay times less 
than the sum of the duration of the excited phase Te 
and two durations of the refractory phase 2Tr, the his- 
tory is non-variable, i.e. there is only one possible profile 
of the two-state process on the interval of time [to — r ; to] , 
where to is the time moment when the latest refractory 
phase has just finished. Therefore, all the waiting times 
are identically distributed and the renewal theory can 
be applied. In this case ip-(£) as function of the delay 
time r is computed straight forwardly and represented 
by Eq. (0) with P_ determined by Eqs. (TO)) , flTTJ) . We 
use the results of the renewal theory [4ll2l| to obtain the 
power spectrum of the stochastic process s(t). The main 
period of the noise-induced oscillations calculated from 
the analytically known power spectrum shows piece-wise 
linear dependence on the delay time. This analytical re- 
sult confirms the similar finding obtained numerically in 
0, [l(| for the FitzHugh-Nagumo system in the excitable 
regime. 

For the delay times larger than (Te + 2Tr) the his- 
tory becomes variable and the distribution density of the 
waiting times is no longer time independent. To han- 
dle renewal processes with history-dependent RTDs the 
equilibrium RTD V'- q (£) Eq. (fT5|) in the sense of the av- 
eraging over all possible histories is introduced. The in- 
tegral Eq. (|17p for ^- q (£) is derived for an arbitrary de- 
lay time t. The solution of this equation is given for 
t e [T E + 2T R ]2T E + 2T R ] by Eq. (2ZJ). The knowl- 
edge of V'i q (0 allows us to calculate the average of any 
given function /(s) of the stochastic process s be means 
of Eq. (f2"5| . This is an exact formula which is valid for 
any renewal process with the dependence on history. Un- 
fortunately, the results of the renewal theory can not be 
used directly to compute e.g. the power spectrum of 
the noise-induced oscillations. However, an analog of the 
mean-field approximation is introduced when the expres- 
sion for the power spectrum from the renewal theory is 
calculated with the equilibrium "0l q (£). 

The analytic results were compared with the results 
of numerical simulation of the bistable system @. It 
is shown that the mean-field power spectrum predicts 
correctly the critical delay time when the first branch 
switching occurs in the dependence of the main period 
vs delay time (Fig. [5Ja),(b)). 

Finally, we demonstrated incoherence maximization 
for positive feedback strength due to delay (Fig. [Ha)). 
Hereby, the degree of incoherence measured by the coeffi- 
cient of variation (|29p is shown to posses a local maximum 
for increasing noise strength D and fixed delay time. For 
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negative feedback strength incoherence maximization is work was supported by EPSRC (UK), 
not observed up to the delay times of (2T E + 2T R ). 
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